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Large eddy simulation^oF 
shock turbulence interaction 


By S. Lee 


1. Motivations and objectives 

The presence of shock waves is an important feature that distinguishes high-speed 
flows from low-speed ones. Understanding the mechanisms of turbulence interacting 
with a shock wave is not only of generic interest, but also of fundamental importance 
in predicting the interactions of turbulent boundary layers with shock waves which 
occur in many practical engineering applications. Many flows of interest have very 
high Reynolds numbers with large ranges of dynamic scales which cannot be cap- 
tured by direct numerical simulation (DNS). Large eddy simulation (LES) is thus 
required where the large scale structures are resolved and the effects of the small, 
unresolved “subgrid” scales on the larger, resolved scales are modeled. The cost of 
resolving the shock wave structure is extremely prohibitive for a strong shock wave 
(Lee et al 1992). Therefore, a numerical method is required which can predict 
changes of flow variables across the shock wave without resolving the structure. 
The method is required to have high order numerical accuracy in smooth regions 
of the flow domain to properly predict the evolution of turbulence. 

A new subgrid-scale (SGS) model was developed by Germano et al. (1991) that 
augments the standard Smagorinsky (1963) eddy viscosity model by replacing the 
ad hoc , flow dependent constant with a coefficient determined by the resolved scales 
in the LES. The coefficient automatically adjusts to the flow conditions. This 
“dynamic” SGS model has been extended by Moin et al. (1991) to the simulation 
of compressible turbulence and passive scalar transport. Recently, Ghosal et al. 
(1992) developed SGS models which eliminate some inconsistencies in the previous 
dynamic models. 

A shock capturing scheme developed by Harten et al. (1987) is essentially non- 
oscillatory (ENO), which can also be constructed to be high order accurate through- 
out the domain. The non-oscillatory nature of the scheme is achieved by choosing a 
interpolating stencil which gives the smoothest evaluation of the derivatives of the 
primitive variables. Shu and Osher (1989) extended the scheme to intei'polate the 
fluxes instead of the primitive variables, which significantly simplifies the scheme. 

The immediate goal of this work is two-folds: to test the performance of LES 
calculation performed with a conservative formulation, a formulation preferred in 
the shock-involved simulations, and to test the performance of the shock capturing 
scheme and validate the scheme against the DNS by Lee et al. (1992) for a weak 
shock wave. 

An immediate extension of the present work is to perform DNS of low Reynolds 
number turbulence interacting with strong shock waves by using the shock capturing 
scheme. Understanding the interaction of high Reynolds number turbulence with 
shock waves through LES is the long range goal of this work. 
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2. Accomplishments 

A nonconservative formulation of the energy equation (solving for internal en- 
ergy) was used to perform large eddy simulations of compressible turbulence in 
Moin et ai (1991) due to its simplicity in implementing SGS models compared to 
the conservative formulation (solving for total energy). In problems with shocks in 
the domain, however, the total energy formulation is preferred due to its conser- 
vative nature. A conservative set of equations for the LES were derived from the 
nonconservative equations derived by Moin et al. (1991). Performance of the con- 
servative formulation was compared with the experiment on decaying grid-generated 
turbulence as well as with the filtered DNS field, which is reported in §2.1. Various 
shock-capturing schemes were tested, and an ENO shock-capturing scheme of Shu 
and Osher (1989) was chosen for the simulation of shock/ turbulence interaction. 
The scheme was tested and validated against the data base generated by DNS of 
weak shock waves, which is reported in §2.2. The results obtained with the ENO 
scheme were within 5% from the DNS results, and used less than 25% of the CPU 
time used in the DNS. 

2.1 LES with a conservative formulation 

In the LES of compressible simulation by Moin et ai (1991), nonconservative 
energy equation was used due to its simple form. The equations used are 
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where <7^, Tik, and qk represent the resolved viscous stress, the subgrid-scale stress, 
and heat flux, 
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qt = p(ttkT - Ukf), 

respectively. The overbar denotes the filtering operation, and the tilde denotes the 
density-weighted filter. The filtering operations always refer to a sharp cutoff (or 
volume- weighted average) in the homogeneous direction of the flow. The subgrid 
scale stress Tik and heat flux qk were determined through the dynamic procedure 
suggested by Germano et al. (1991). Derivation of the energy equation for the LES 
using the conservative formulation introduces terms which are complex to model, 
for example, the subgrid scale convection of the total energy (sum of thermal and 
kinetic energy) needs to be modeled. ^ 

However, if the resolved total energy (Et) is defined as 
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FIGURE 1. Evolution of resolved turbulent kinetic energy. • filtered experi- 
ment (Comte-Bellot & Corrsin 1971), LES with conservative formulation, 

LES with nonconservative formulation. 


Et = p(c v T + ukUk/2) 

rather than p(c v T + ti^fe/2), the conservative energy equation can be derived from 
the nonconservative equations presented above, and the extra complexities can be 
circumvented. The resulting energy equation becomes 


dEx dEx^k 
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The terms to be modeled remain the same and are parametrized as 


r i ;=^C s A 2 |5|5* fcl 
tf 2 = ^C/A 2 |S| 2 , 
qt =pDA 2 \S\T k , 

where q 2 = r mm and r* k = r lk — q 2 S lk / 3. The detailed procedures of dynamically 
determining coefficients, Cs, Cj , and D, are described in Moin et al (1991). The 
concept of least squares suggested by Lilly (1992) replaced the rather arbitrary 
procedure of determining scalar coefficients, Cs and Z), from the tensor or vector 
relations. 

The experiment on the decay of grid-generated turbulence with Re\ = 71.6 
(Comte-Bellot and Corrsin 1971) was simulated as a temporal decay with 32x32x32 
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Figure 2. Evolution of resolved density and temperature fluctuations. Conser- 
vative formulation: density, temperature, nonconservative formulation: 

• density, * temperature. 

grid points with different formulations. Decay of the resolved scale TKE is com- 
pared in Figure 1. The evolution is insensitive to the choice of the formulation. 
Evolutions of rms density and temperature fluctuations are compared in Figure 2. 
(The fluctuations in thermodynamic quantities were initially set to zero.) Even 
though the evolution of velocity field does not depend significantly on the choice of 
the formulation, the evolution of density fluctuation depends on this choice. 

To check if this dependence is due to the initial transient of the fluctuation density 
field, LES’s of decaying compressible turbulence (Re\ = 35, M t = 0.61, where eddy- 
shocklets were found in DNS) were performed. Initial conditions of the simulation 
were taken as the Fourier filtered DNS field onto 32 x 32 x 32 grids at t/r* ~ 0.85, 
when turbulence is fully developed in every respect. Figure 3 shows evolution of 
rms density fluctuation from the two formulations, which shows good agreement 
between the two and with the filtered DNS. 

The reason for the different evolution of density field shown in Figure 2 is not 
yet fully understood. Some possible reasons are: (1) aliased evaluation of temper- 
ature in the conservative formulation is subtracted from £r), (2) different 

initial transient to set up temperature fluctuations. Before conducting the LES 
of shock/turbulence interaction, validity of the formulation and the SGS models 
should be established. 

2.2 Development of an ENO scheme 

For a strong shock wave, resolving the shock structure by solving Navicr-Stokes 
equations is irrelevant because Navier- Stokes equations are no longer valid inside 
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FIGURE 3. Evolution of resolved density fluctuation. • filtered DNS (Lee et al. 

1992), LES with conservative formulation, — -- LES with nonconservative 

formulation. 


the shock wave with M\ > 2.0. (Sherman 1955) Therefore, there is no point to 
resolving the structure of the strong shock wave with the expense of huge computa- 
tional cost. A shock-capturing scheme is implemented for the simulation of strong 
shock wave/ turbulence interaction. In simulations of shock/turbulence interaction, 
a shock capturing scheme is required to be high order accurate throughout the 
computational domain (to properly simulate the turbulence evolution) as well as to 
be able to give smooth and accurate transition across the shock wave. The shock- 
capturing scheme of the choice is the essentially non- oscillatory (ENO) scheme, 
which can be constructed up to any order of accuracy. The ENO scheme used in 
this work is based on the Lax- Friedrichs scheme with interpolation of fluxes (Shu 
and Osher, 1989). Some modifications are made to the basic scheme to improve the 
solution accuracy and to enhance the code performance. 

For the sake of completeness, the basic ENO scheme is briefly described in the fol- 
lowing. The three-dimensional, compressible Navier-Stokes equation can be written 
as 


q* + f(q)x + g(q)» + h(q)* = v. 


where, 
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q = ( p,pui,pu 2 ,pu 3 ,E T ) T , 
f(q) = «iq + p( o, l, 0, 0, «j ) T , 
g(q) = u 2 q +p(°> 0 * 1 > 0 > u 2 ) T i 
h(q) = « 3 q + pC0,0,0,l,u 3 ) T , 

and V is the molecular diffusion terms. 

The ENO adaptive-stencil procedure is applied only to the convection terms (the 
left-hand side of the above equation). The diffusion terms are approximated by 
the sixth-order Pade scheme, which is the base scheme used in the shock- resolving 
direct numerical simulation (Lee et al. 1992). The time advancement is performed 
by a third order Runge-Kutta method. Three convection terms are approximated 
dimension- by-dimension: for example, when approximating f(q) x , y and z are fixed. 
The core algorithm is then a one-dimensional ENO approximation to f(q) c . 

The flux vector f(q) x can be approximated either by component by component 
or in characteristic directions, but for simplicity of implementation, the compo- 
nent treatment was chosen. To obtain nonlinear stability by upwinding, f(q) r is 
decomposed into 


f(q)« = f + (q)x + f~(q)x 

with 

^(q)* = (f(q)x ± aq)/2, 

where a = max(|ui| + c) (c is the sound speed) is the largest eigenvalue in absolute 
value of the Jacobian df/dq along the relevant z-line. The decomposition guaran- 
tees that df^/dq has positive/negative eigenvalues. Therefore, upwinding (to be 
discussed) is the same for all components of the flux. 

Suppose / is a component of the flux f ± (q). /(q)* I s written as a conserved flux 

difference 


/■I-,, - £- x v>*i - fi-O’ 

where the numerical flux fj+ 1/2 approximates /i(xj +1 / 2 ) to a high order with h(x) 
defined by 


/(«(*)) = 


Ax 



In order to evaluate the numerical flux /j+ 1 / 2 , h(x) need not to be explicitly con- 
structed. The numerical fluxes can be expressed in terms of the undivided differ- 
ences of f(u(x)) as 


r 

fj+i/2 = c (* 
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with i being the left-most point in the stencil used to approximate /,+ x / 2 , and 
c(s,m) and the undivided difference /[?, A:] being defined by 


and 


c(s, m) 


1 

(m + 1)! 


3+m i+m 

E IB-*) 

*=» iz: 


f[h o] = /(«>), 

f[j,k] = f[] + l,k - 1] - f{j,k - 1], k = 1 


respectively. 

The adaptive stencil (represented by i ) is determined by the following procedure. 

(1) Start with i = j or i = j + 1 according to the local upwinding direction. 

(2) For k~ th level of approximation, 

i = i — 1 if |/[i, £]| > |/[z — 1, fc]| for k = 1, ...,r. 

The ENO scheme was implemented in the code developed for the simulation of 
shock/ turbulence interaction. Simulation of spatially decaying turbulence was per- 
formed (1) to validate the scheme and (2) to determine the required order of the 
ENO scheme to reproduce the quality of the corresponding DNS. In the validation 
procedure, a totally unexpected accuracy degeneracy was noticed with the increase 
of the ENO scheme’s order. As the order of the ENO increases, the solution devi- 
ates further from the DNS solution. Similar accuracy degeneracy was reported by 
Rogerson and Meiburg (1990) while they were refining the grid. This degeneracy 
was due to the choice of linearly unstable stencils during the adaptive ENO proce- 
dure (Shu 1990). This phenomenon is more pronounced for the higher order scheme 
with more refined grid. A biasing toward the central difference stencil (Shu 1990) 
was implemented to avoid the accuracy degeneracy in smooth regions. With this fix, 
by raising the order of accuracy, more accurate results could be obtained. Results 
comparable to the corresponding DNS (with sixth-order Pade scheme) was obtained 
when the accuracy of the ENO scheme was sixth order. Beyond the sixth-order, 
the effect of the increased order was found negligible in this case. 

The operation count and, accordingly, the CPU time are increased a few times 
compared to the standard DNS code when the ENO scheme is implemented through- 
out the domain. The region where the ENO scheme plays a significant role is quite 
localized, usually less than a tenth of the computational domain in DNS of weak 
shock waves, and is expected to be even smaller for stronger shock waves. The 
idea of applying the ENO scheme only where it plays an active role and switch- 
ing elsewhere into the usual Pade central difference was tested in pursuit of saving 
unnecessary computer usage. The local region of ENO application can even be 
specified a priori to be a zone around the shock since the simulation is performed 
in a coordinate system fixed on the mean shock wave. The concept of local applica- 
tion of the ENO scheme was validated against fully resolved DNS results both for 
the case of two-dimensional shock/turbulence interaction with Mi = 1.2, M t = 0.07 
and the case of three-dimensional spatially decaying turbulence with Mt = 0.51 and 
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FIGURE 4. Evolution of velocity fluctuations in three-dimensional shock/turbulence 
interaction. Lines denote DNS results and symbols denote results using the ENO 
scheme, uj ( , • ), w 2 ( > x )> u 3 ( , + )• 

Re\ = 25. For the two-dimensional shock/turbulence case with M\ = 1.2, CPU 
time required was reduced to a third of the corresponding DNS and is expected to 
be larger for the stronger shock waves. 

The next validation was to simulate the three-dimensional shock/ turbulence in- 
teraction using the ENO scheme near the shock and the Pade scheme elsewhere. 
Significant additional decay of TKE was found in the shock vicinity where the ENO 
scheme is applied, which was not pronounced in two-dimensional shock/turbulence 
simulations. The excessive TKE dissipation in the ENO zone was suspected to be 
from the additional dimension (x 3 -direction) for the ENO procedure. In order to 
investigate the effect of the dimensions, the ENO scheme was applied only in the 
shock-normal direction for the two-dimensional shock/turbulence simulation. More 
accurate results were obtained with less CPU time used. For M\ = 1.2, the ENO 
scheme required only a quarter of the CPU time used in the corresponding DNS. 

The final ENO scheme developed in this work interpolates (up to sixth-order) 
fluxes based on Lax- Friedrichs method with biasing the interpolation stencil toward 
the central differencing. The ENO procedure is applied only to the convection terms 
in the shock-normal direction. The ENO scheme is made only active near the shock 
wave and switched into the sixth-order Pade scheme elsewhere. 

The final validations were made for two cases: (1) three-dimensional turbulence 
interacting with a weak shock wave {M\ = 1.2) and (2) two-dimensional turbulence 
interacting with a relatively strong shock wave ( M\ — 2.0). The first case was 
chosen to check the performance of the ENO scheme in three-dimensional case, and 
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FIGURE 5. Evolution of vorticity fluctuations in three-dimensional shock/turbulence 
interaction. Lines denote DNS results and symbols denote results using the ENO 
scheme, u \ ( , • ), u; 2 ( , * )> ( , 4- ). 

the result was compared with the corresponding DNS result (Lee et al. 1992). The 
second case was chosen to check the accuracy of the ENO scheme for strong shock 
waves. For the purpose of comparison, a shock-resolving DNS of two-dimensional 
shock/turbulence interaction was also performed. 

As a validation of the ENO scheme for three-dimensional shock/turbulence in- 
teraction, reproduction of the DNS results of Case C in Lee et al. (1992) was tried 
where M\ = 1.2, M t = 0.095, and Re\ = 12 just upstream of the shock. In the 
DNS, 193 x 64 x 64 grids (with significant grid refinement near the shock) and 160 
Cray/Y-MP CPU hours were used to give marginally converged statistics. In the 
ENO scheme, 129 x 64 x 64 uniform grids and 40 Cray/Y-MP CPU hours were used 
for the same quality of statistics. 

Figure 4 compares the evolution of velocity fluctuations obtained by DNS and 
ENO in three-dimensional shock/turbulence interaction. Except for the position of 
the mean shock wave, the amplification and free evolution of velocity from DNS 
are well reproduced by the ENO scheme. Since the ENO scheme is activated in 
the shocked zone, the peak rms fluctuations inside the shock wave due to the shock 
intermittency are not expected to agree with the DNS. 

Figure 5 compares vorticity component evolution obtained by DNS and ENO. 
Once the shift of the mean shock position is considered, predictions by the ENO 
scheme are in excellent agreement with the DNS results. Other turbulence quan- 
tities, such as fluctuations in dilatation, density, and temperature, were in good 
agreement with the DNS results outside the zone occupied by the shock wave. 
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FIGURE 6(a) . Evolution of velocity fluctuations in two-dimensional shock/ turbulence 
interaction. Lines denote DNS results and symbols denote results using the ENO 
scheme, uj ( , • ), U2 ( , x )* 
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Figure 6(b) . Evolution of vorticity fluctuations in two-dimensional shock/ turbulence 
interaction. DNS, ENO. 
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As a validation of the ENO scheme for shock/turbulence interaction at high Mach 
numbers, results of the ENO scheme on two-dimensional shock/turbulence interac- 
tion was compared with the corresponding DNS, where Mj = 2.0 and M t = 0.07 at 
the inflow. The minimum grid spacing required by the ENO scheme was five times 
larger than that used in the DNS. Figure 6(a) and (b) compare velocity and vor- 
ticity fluctuations from the ENO scheme with those from the DNS. The jumps and 
evolutions behind the shock wave were reproduced by the ENO scheme with much 
less cost than the DNS. The results from the ENO scheme are indeed dependent on 
the resolution near the shock wave, and a grid-refinement test was performed before 
comparing the results. In order to accurately predict the interaction of turbulence 
with the stronger shock waves, where no DNS results are available, a grid in the 
shock vicinity should be refined until the result is no longer dependent on further 
grid refinement. 

3. Future plans 

Further tests of the LES with the conservative formulation against the filtered 
DNS and LES with the nonconservative formulation will continue. Extension of 
Ghosal et al.'s (1992) model (a constrained integral equation approach and a SGS 
model using transport equation for SGS turbulent kinetic energy) to compressible 
turbulence will also be performed. 

In most DNS and LES of compressible turbulence, aliasing errors were controlled 
rather than clearly removed. In Lee et al. (1992), it was proposed to solve a 
continuity equation in terms of specific volume instead of density to completely 
remove aliasing errors. A simulation with this formulation will be executed. The 
simulation can be used as a reference to investigate the accuracy of the existing 
formulations. 

Interaction of three-dimensional isotropic turbulence with strong shock waves 
will be performed to study the effect of the shock strength. Once SGS models and 
formulations are developed and validated, LES of turbulence interacting with shock 
waves of various strengths will also be performed. 
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